
########################## Start: Copied from the main regressions ##################
################################### Packages ##########################################
library(tidyverse)
library(fixest)
library(stargazer)

################################### Set Working Directory ##########################################

# Getting the path of your current file
setwd(dirname(rstudioapi::getSourceEditorContext()$path))

################################### functions ##########################################
# IHS Transformation
ihs=function(x) log(x+sqrt(1+x^2)) # to revers ihs use sinh
# Demeans and IHS Transforms
idm=function(x) ihs(x)-mean(ihs(x),na.rm=T)
# Demeans
dm=function(x) x-mean(x,na.rm=T)



################################### Data ##########################################
### main data
up = read.csv("../Data/data_main_final.csv",stringsAsFactors = F)%>%
  ### The following codes IHS transforms the variables and demeans them at the stock level. The demeaning is only relevant for the interacted variables
  group_by(stock_id)%>% # for demeaning at the stock level
  arrange(year)%>% # for calculating the lagged variables
  mutate(
    ffmsy = ihs(ffmsy), # Use only IHS because of the simulation (below). The demeaning is only important for the interacted variables, for the non-iteracted variables nothing changes
    gdp=idm(gdp_const),
    gdp2=idm(gdp_const)^2,
    tfp=idm(rtfpna),
    fao_ag_credit = idm(fao_agriculture_credit),
    fao_to_credit = idm(fao_total_credit),
    credit=idm(domestic_credit),
    credit_instrument = idm(credit_instrument),
    interest_factor=dm(1/(1 + lending_interest/100)),
    interest_factor = ifelse(inflation_cp<=30,interest_factor,NA), # This removes observations during hyper inflation
    inflation_factor=dm(1/(1+inflation_cp/100)),
    inflation_factor=ifelse(inflation_cp<=30,inflation_factor,NA), # This removes observations during hyper inflation
    bbmsy1 = idm(lag(bbmsy,1))
  )%>%
  as.data.frame()


### ram data
ram = read.csv("../Data/data_ram_final.csv",stringsAsFactors = F)%>%
  ### The following codes IHS transforms the variables and demeans them at the stock level. The demeaning is only relevant for the interacted variables
  group_by(country, stock_id, scientific_name)%>%
  arrange(year)%>% # for calculating the lagged variables
  mutate(ffmsy = ihs(ffmsy), # Fishing mortality
         gdp=idm(gdp_const), # GDP (World Bank)
         gdp2=idm(gdp_const)^2, 
         tfp = idm(rtfpna), # Total factor productivity (Penn)
         fao_ag_credit = idm(fao_agriculture_credit), # Credit to agriculture, fisheries, forestry from the FAO
         fao_to_credit = idm(fao_total_credit), # Total credit from the FAO
         credit=idm(domestic_credit), # Domestic credit to private sector (% of GDP)
         interest=idm(lending_interest), # Lending interest IMF
         inflation=idm(inflation_cp), # Inflation from World Bank
         interest_factor=dm(1/(1 + lending_interest/100)),
         inflation_factor=dm(1/(1+inflation_cp/100)),
         bbmsy1 = idm(lag(bbmsy,1)), # Fish stocks 
         species_cat_id =as.numeric(as.factor(isscaap_group)) # Species group
  )%>%
  as.data.frame()

########################## End: Copied from the main regressions ##################

############## Results for Appendix A.12 Property rights and the Rule of Law #########################
### Defines the new property rights variable ###############
up_rol = up%>%
  mutate(pr = pr*wbgi_rle_stine)


### 1. Baseline without interaction term
c1 = feols(ffmsy ~ credit + pr 
           + gdp + gdp2 
           | stock_id + country + year,
           cluster = ~ country + year + scientific_name,
           data = up_rol)

summary(c1)



### 2. Main
c2 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up_rol)

summary(c2)



### 3. trend
c3 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year + country[year],
           cluster = ~ country + year + scientific_name, 
           data = up_rol)

summary(c3)


#### 4. IV (the instrument starts in 1977)

c4 = feols(ffmsy ~  pr 
           + gdp + gdp2 
           | stock_id + country + year
           |credit + credit:pr ~ credit_instrument + credit_instrument:pr ,
           cluster = ~ country + year + scientific_name, 
           data = up_rol) # Note that the credit instrument data are only available after 1977, so many observations drop out


summary(c4,stage = c(1,2))
fitstat(c4,type = "ivwald")



#### 5. OECD
# These are all OECD countries but note that the land-locked countries are not in the data (Austria, Czech Republic, Hungary,...)
oecd = c("Australia" , "Austria" , "Belgium" , "Canada" , "Chile" , "Colombia" , 
         "Costa Rica" , "Czech Republic" , "Denmark" , "Estonia" , "Finland" , "France" , "Germany" , 
         "Greece" , "Hungary" , "Iceland" , "Ireland" , "Israel" , "Italy" , "Japan"  ,
         "Latvia" , "Lithuania" , "Luxembourg" , "Mexico" , "Netherlands" , "New Zealand" , "Norway" ,
         "Poland" , "Portugal" , "Slovakia" , "Slovenia" , "South Korea" , "Spain" , "Sweden" , "Switzerland" , "Turkey" , 
         "United Kingdom" , "United States")



c5 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up_rol%>%filter(country%in%oecd))

summary(c5)


##### 6. Balanced
years = 1960:2012
# number of years
ny = length(years)

# credit
c6 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up_rol%>% # This code restricts the stocks to those with complete data. 
             filter(year%in%years)%>%
             group_by(stock_id,country)%>%
             filter(n()==ny) # this removes all stocks with less obsevrations than the number of years. Note that all NA values for the ffmsy variable were removed earlier
)

summary(c6)

##### 7. Top10
### Orders species by total catch value between 2002 and 2012. We restrict this to this period to reduce the influence of random fluctuations (less years) and the influence of time series length (more years)
top = up_rol%>%
  filter(year%in%2002:2012)%>%
  group_by(scientific_name)%>%
  summarise(value = sum(catch*price))%>% # Sums the catch values over the specified time period
  arrange(desc(value))


### selects top 10% of species
# selects all species names
spe = top$scientific_name
# number of species in the top 10 % (rounding is necessary to reduce the number to integers)
p10 = round(length(spe)*0.1,0) 
# top 10 % species
top10 = spe[1:p10]

# regression
c7 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year,
           cluster = ~ country + year + scientific_name, 
           data = up_rol%>%filter(scientific_name%in%top10))

summary(c7)




#### 8. Ram
c8 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = ram%>%
             mutate(pr = pr*wbgi_rle_stine))

summary(c8)


#### 9. Ag. Credit

c9 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 + fao_to_credit 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up_rol%>%mutate(credit = fao_ag_credit))

summary(c9)




###### print results
etable(c1,c2,c3,c4,c5,c6,c7,c8,c9,
       digits = 3,
       digits.stats = 3,
       drop = c("gdp","gdp2","unemployment_nat","rtfpna","fao_to_credit"),
       dict=c("pr" ="Property rights",
              "credit" = "Credit",
              "stock_id" = "Stock FE",
              "country" = "Country FE",
              "year" ="Year FE",
              "credit x pr" = "Credit $\\times$ Property rights",
              "credit:pr" = "Credit $\\times$ Property rights",
              "Year FE (Country FE)" = "Country trend" ),
       replace =T,
       order = c("!times","Property rights","Credit"),
       headers = c("Baseline","Main","Trend","IV","OECD","Balanced","Top10","RAM","Ag. Credit"),
       depvar = F,
       drop.section = c("slopes"),
       fixef_sizes =T,
       file = "../Results/results-credit-rol_final.tex")

############## Appendix A.13 Total factor productivity and the resource stock #########################

### 1. baseline
c1 = feols(ffmsy ~ credit + pr 
           + gdp + gdp2 + bbmsy1 + tfp
           | stock_id + country + year,
           cluster = ~ country + year + scientific_name,
           data = up)

summary(c1)



### 2. main

c2 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 + bbmsy1 + tfp
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up)

summary(c2)



### 3. trend
c3 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 + bbmsy1 + tfp
           | stock_id + country + year + country[year],
           cluster = ~ country + year + scientific_name, data = up)

summary(c3)


#### 4. IV (the instrument starts in 1977)
c4fs = feols(ffmsy ~  pr 
             + gdp + gdp2 + bbmsy1 + tfp
             | stock_id + country + year
             |credit  ~ credit_instrument ,
             cluster = ~ country + year + scientific_name, 
             data = up%>%
               filter(year>=1977))

summary(c4fs, stage = 1)
fitstat(c4fs,type = "ivwald")


c4 = feols(ffmsy ~  pr 
           + gdp + gdp2 + bbmsy1 + tfp
           | stock_id + country + year
           |credit + credit:pr ~ credit_instrument + credit_instrument:pr ,
           cluster = ~ country + year + scientific_name, 
           data = up%>%
             filter(year>=1977))


summary(c4,stage = c(1,2))
fitstat(c4,type = "ivwald")



#### 5. OECD
oecd = c("Australia" , "Austria" , "Belgium" , "Canada" , "Chile" , "Colombia" , 
         "Costa Rica" , "Czech Republic" , "Denmark" , "Estonia" , "Finland" , "France" , "Germany" , 
         "Greece" , "Hungary" , "Iceland" , "Ireland" , "Israel" , "Italy" , "Japan"  ,
         "Latvia" , "Lithuania" , "Luxembourg" , "Mexico" , "Netherlands" , "New Zealand" , "Norway" ,
         "Poland" , "Portugal" , "Slovakia" , "Slovenia" , "South Korea" , "Spain" , "Sweden" , "Switzerland" , "Turkey" , 
         "United Kingdom" , "United States")

# credit
c5 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 + bbmsy1 + tfp
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, data = up%>%filter(country%in%oecd))

summary(c5)


##### balanced
years = 1960:2012
ny = length(years)

# credit
c6 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 + bbmsy1 + tfp
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up%>%
             filter(year%in%years)%>%
             group_by(stock_id,country)%>%
             filter(n()==ny))

summary(c6)


### top species in terms of catch value between 2002 and 2012
top = up%>%
  filter(year%in%2002:2012)%>%
  group_by(scientific_name)%>%
  summarise(value = sum(catch*price))%>% 
  arrange(desc(value))


# top 10%
spe = top$scientific_name
p10 = round(length(spe)*0.1,0)
top10 = spe[1:p10]


c7 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 + bbmsy1 + tfp
           | stock_id + country + year,
           cluster = ~ country + year + scientific_name, 
           data = up%>%
             filter(scientific_name%in%top10))

summary(c7)




#### 8. ram
c8 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 + bbmsy1 + tfp
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, data = ram)

summary(c8)


#### 9. fao credit

c9 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 + bbmsy1 + tfp+ fao_to_credit 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up%>%mutate(credit = fao_ag_credit))

summary(c9)

###### print results
etable(c1,c2,c3,c4,c5,c6,c7,c8,c9,
       digits = 3,
       digits.stats = 3,
       drop = c("gdp","gdp2","unemployment_nat","rtfpna","fao_to_credit","bbmsy1","tfp"),
       dict=c("pr" ="Property rights",
              "credit" = "Credit",
              "stock_id" = "Stock FE",
              "country" = "Country FE",
              "year" ="Year FE",
              "credit x pr" = "Credit $\\times$ Property rights",
              "credit:pr" = "Credit $\\times$ Property rights",
              "Year FE (Country FE)" = "Country trend" ),
       replace =T,
       order = c("!times","Property rights","Credit"),
       headers = c("Baseline","Main","Trend","IV","OECD","Balanced","Top10","RAM","Ag. Credit"),
       depvar = F,
       drop.section = c("slopes"),
       fixef_sizes =T,
       file = "../Results/results-credit-tfp_final.tex")


########################## Appendix A.14 Excluding multinational stocks #############################################################
### Robustness excluding multinational stocks
upm = up%>%
  filter(multinational==0) 



### 1. Baseline without interaction term
c1m = feols(ffmsy ~ credit + pr 
           + gdp + gdp2 
           | stock_id + country + year,
           cluster = ~ country + year + scientific_name,
           data = upm)

summary(c1m)



### 2. Main
c2m = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = upm)

summary(c2m)



### 3. trend
c3m = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year + country[year],
           cluster = ~ country + year + scientific_name, 
           data = upm)

summary(c3m)


#### 4. IV (the instrument starts in 1977)

c4m = feols(ffmsy ~  pr 
           + gdp + gdp2 
           | stock_id + country + year
           |credit + credit:pr ~ credit_instrument + credit_instrument:pr ,
           cluster = ~ country + year + scientific_name, 
           data = upm) # Note that the credit instrument data are only available after 1977, so many observations drop out


summary(c4m,stage = c(1,2))
fitstat(c4m,type = "ivwald")



#### 5. OECD
# These are all OECD countries but note that the land-locked countries are not in the data (Austria, Czech Republic, Hungary,...)
oecd = c("Australia" , "Austria" , "Belgium" , "Canada" , "Chile" , "Colombia" , 
         "Costa Rica" , "Czech Republic" , "Denmark" , "Estonia" , "Finland" , "France" , "Germany" , 
         "Greece" , "Hungary" , "Iceland" , "Ireland" , "Israel" , "Italy" , "Japan"  ,
         "Latvia" , "Lithuania" , "Luxembourg" , "Mexico" , "Netherlands" , "New Zealand" , "Norway" ,
         "Poland" , "Portugal" , "Slovakia" , "Slovenia" , "South Korea" , "Spain" , "Sweden" , "Switzerland" , "Turkey" , 
         "United Kingdom" , "United States")



c5m = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = upm%>%filter(country%in%oecd))

summary(c5m)


##### 6. Balanced
years = 1960:2012
# number of years
ny = length(years)

# credit
c6m = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = upm%>% # This code restricts the stocks to those with complete data. 
             filter(year%in%years)%>%
             group_by(stock_id,country)%>%
             filter(n()==ny) # this removes all stocks with less obsevrations than the number of years. Note that all NA values for the ffmsy variable were removed earlier
)

summary(c6m)

##### 7. Top10
### Orders species by total catch value between 2002 and 2012. We restrict this to this period to reduce the influence of random fluctuations (less years) and the influence of time series length (more years)
top = upm%>%
  filter(year%in%2002:2012)%>%
  group_by(scientific_name)%>%
  summarise(value = sum(catch*price))%>% # Sums the catch values over the specified time period
  arrange(desc(value))


### selects top 10% of species
# selects all species names
spe = top$scientific_name
# number of species in the top 10 % (rounding is necessary to reduce the number to integers)
p10 = round(length(spe)*0.1,0) 
# top 10 % species
top10 = spe[1:p10]

# regression
c7m = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year,
           cluster = ~ country + year + scientific_name, 
           data = upm%>%filter(scientific_name%in%top10))

summary(c7m)





###### print results
etable(c1m,c2m,c3m,c4m,c5m,c6m,c7m,
       digits = 3,
       digits.stats = 3,
       drop = c("gdp","gdp2","unemployment_nat","rtfpna","fao_to_credit"),
       dict=c("pr" ="Property rights",
              "credit" = "Credit",
              "stock_id" = "Stock FE",
              "country" = "Country FE",
              "year" ="Year FE",
              "credit x pr" = "Credit $\\times$ Property rights",
              "credit:pr" = "Credit $\\times$ Property rights",
              "Year FE (Country FE)" = "Country trend" ),
       replace =T,
       order = c("!times","Property rights","Credit"),
       headers = c("Baseline","Main","Trend","IV","OECD","Balanced","Top10"),
       depvar = F,
       drop.section = c("slopes"),
       fixef_sizes =T,
       file = "../Results/main_final_multinational.tex")







############## Appendix A.15 Discount Factor #########################

### 1. baseline
i1 = feols(ffmsy ~ interest_factor + pr 
           + gdp + gdp2 + inflation_factor
           | stock_id + country + year,
           cluster = ~ country + year + scientific_name,
           data = up)

summary(i1)



### 2. main

i2 = feols(ffmsy ~ interest_factor + pr + interest_factor:pr 
           + gdp + gdp2 + inflation_factor
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up)

summary(i2)



### 3. trend
i3 = feols(ffmsy ~ interest_factor + pr + interest_factor:pr 
           + gdp + gdp2 #+ inflation_factor
           | stock_id + country + year + country[year],
           cluster = ~ country + year + scientific_name, data = up)

summary(i3)




##### balanced
years = 1960:2012
ny = length(years)

# interest_factor
i4 = feols(ffmsy ~ interest_factor + pr + interest_factor:pr 
           + gdp + gdp2 + inflation_factor
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up%>%
             filter(year%in%years)%>%
             group_by(stock_id,country)%>%
             filter(n()==ny))

summary(i4)





#### 8. ram
i5 = feols(ffmsy ~ interest_factor + pr + interest_factor:pr 
           + gdp + gdp2 + inflation_factor
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, data = ram)

summary(i4)



###### print results
etable(i1,i2,i3,i4,i5,
       digits = 3,
       digits.stats = 3,
       drop = c("gdp","gdp2","unemployment_nat","rtfpna","fao_to_credit","inflation_factor"),
       dict=c("pr" ="Property rights",
              "interest_factor" = "Discount factor",
              "stock_id" = "Stock FE",
              "country" = "Country FE",
              "year" ="Year FE",
              "Year FE (Country FE)" = "Country trend" ),
       replace =T,
       order = c("!times","Property rights","Credit"),
       headers = c("Baseline","Main","Trend","Balanced","RAM"),
       depvar = F,
       drop.section = c("slopes"),
       fixef_sizes =T,
       file = "../Results/results-interest_final.tex")

